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ABSTRACT 


A  numerical  model  for  heat  transfer  in  laminar  duct  flows  has  been  developed 
using  the  finite  difference  method  to  explore  the  significance  and  extent  of  “back- 
conduction”  at  low  Peclet  numbers.  The  calculations  have  been  carried  out  for  flows 
between  parallel  plates  and  in  circular  tubes  by  using  different  Peclet  numbers  in  the 
range  of  0.05  to  100.  For  both  situations  constant  heat  flux  and  constant  wall  temperature 
boimdary  conditions  were  used.  The  validity  of  the  results  has  been  checked  by 
comparison  with  some  existing  results  in  the  literature,  and  extended  to  a  wider  range  of 
parameters  including  conjugate  wall  conduction  effects.  The  results  are  presented  for  bulk 
mean  temperature  variation,  Nusselt  number  behavior,  and  energy  absorbed  before  the 
heated  section,  for  cases  with  and  without  wall  conduction.  Such  axial  conduction  effects 
may  be  an  important  feature  in  the  thermal  characterization  of  microtubes,  which  are  to 
be  used  in  microheat  exchangers. 
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1.  INTRODUCTION 


The  heat  transfer  problem  of  laminar  fluid  flow  in  ducts  which  is  known  as  the 
Graetz  problem  has  many  applications  in  technology  and  has  been  studied  extensively 
since  Graetz  (1885).  The  classical  Graetz  problem  considers  the  forced  convection  heat 
transfer  of  the  fluids  flowing  in  ducts  neglecting  axial  conduction  effects  in  the  fluid  or 
the  wall. 

This  study  is  on  the  axial  conduction  effects  of  the  flows  whose  Peclet  numbers 
range  is  from  those  of  liquid  metals  where  axial  conduction  may  not  be  neglected,  to 
those  of  oils  where  the  axial  conduction  has  almost  no  effect  on  the  temperature 
distribution  inside  the  ducts.  Developing  micromachines  nowadays  and  in  the  future,  the 
axial  conduction  effects  may  be  a  very  important  feature  in  microheat  exchangers. 

The  effects  of  axial  conduction  on  the  Hagen-Poiseuille  flows  between  parallel  plates  and 
in  circular  cylinders  for  constant  wall  temperature  and  constant  heat  flux  have  been 
studied  in  this  thesis.  The  governing  equations  have  been  solved  using  the  finite 
difference  numerical  scheme.  The  validity  of  the  results  has  been  checked  by  comparison 
with  some  existing  results  in  the  literature  and  extended  to  a  wider  range  of  parameters 
including  conjugate  wall  conduction  effects.  The  results  are  presented  for  bulk  mean 
temperature  variation,  Nusselt  number  behavior,  and  energy  absorbed  before  the  heated 
section,  for  cases  with  and  without  wall  conduction. 


1 
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n.  BACKGROUND 


A.  HEAT  TRANSFER  IN  LAMINAR  DUCT  FLOWS 

Consider  a  flow  in  a  circular  tube  where  fluid  enters  the  tube  with  a  uniform 
velocity.  The  viscous  effects  are  important  in  this  flow  and  a  boundary  layer  develops  as 
X  increases.  The  flow  is  “fully  developed”  in  the  region  where  the  velocity  gradient 
doesn’t  change  anymore  with  increasing  x  (Figure  1). 


Figure  1.  Laminar  fully  developed  flow 

The  Reynolds  Number  for  the  duct  flows  is  delSned  as 

(2.1) 

y 

where  V  is  the  mean  fluid  velocity,  v  is  kinematic  viscosity,  v  =  — ,  and  is  the 

P 

hydraulic  diameter,  ^  ^  Cross _Sectionai _Area .  For  the  laminar  duct  flows, 

^  Wetted  _Perimeter 

Rej^  <  2300 . 

The  Prandtl  number  relates  the  temperature  distribution  to  the  velocity  distribution 
and  defined  as 

Pr  =  -  (2.2) 

a 
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where  a  is  the  thermal  dififusivity  .  If  the  Prandtl  number  is  one,  the  velocity  and  the 
temperature  profiles  develop  together  and  at  the  same  rate.  The  Prandtl  number  is: 

Pr  « 1  for  the  liquid  metals 

Pr  =  1  for  the  gases 

Pr»l  for  the  oils 

For  the  liquid  metals,  where  the  Pr«l,  the  energy  diffusion  rate  is  much  more 
than  the  momentum  diffusion  rate.  It  is  opposite  for  the  oils,  in  which  Pr»l  and  the 
velocity  profile  develops  faster  than  the  temperature  profile  in  this  case.  The  value  of  the 
Prandtl  number  strongly  affects  the  relative  growth  of  the  velocity  and  thermal  boundary 
layers. 

Peclet  number  is  defined  as 

Pe=Re.Pr  (2.3) 

and  it  is  a  measure  of  the  quantity  of  the  axial  heat  conduction  effects  in  the  fluid.  The 
axial  conduction  is  assumed  negligible  for  Pe>10  and  the  axial  conduction  term  can  be 
assumed  small  in  the  governing  equation  for  this  case.  But  the  axial  conduction  effects 
can  be  significant  when  the  Peclet  number  is  smaller.  The  purpose  of  this  numerical  study 
is  to  show  the  importance  of  the  axial  conduction  effects  and  the  heat  absorption  upto  the 
entrance  firom  where  the  heating  starts  when  the  Peclet  number  is  small. 

The  Nusselt  number  is  defined  as 

Nu  =  ^  (2.4) 

k  . 

where  k  is  the  thermal  conductivity  of  the  fluid.  It  is  a  non-dimensional  number  where  the 
convection  coefficient  h  is  calculated. 

The  laminar  duct  flow  being  considered  is  assumed  to  be  hydrodynamically  fully 
developed  before  any  heating  effects  are  considered.  The  velocity  profile  in  the  laminar 
fully  developed  region  is: 

For  the  circular  tube: 


u(r)  =  2V 


(2.5) 
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For  the  parallel  plates: 


u(y)  =  -v| 


1- 


(2.6) 


where  r  and  y  are  measiired  from  the  centerline,  ro  is  tube  radius  and  h  is  the  distance 
from  the  centerline  to  the  plates.  [Ref.l] 

1.  The  Heat  Transfer 


The  heat  transfer  for  a  duct  flow  can  be  expressed  using  Newton’s  law  of  cooling, 
q  =h(Tw-Tni),  where  h  is  the  convection  heat  transfer  coefficient,  Tw  is  the  wall 
temperature  and  Tm  is  the  mean  temperature  of  the  fluid,  where  Tm  is 


Jpuc^TdA 

- .  (2.7) 

me,. 

For  the  constant  Cv  and  incompressible  flow  through  the  circular  tube,  Tm  is 


T„=:A]^Trdr  (2.8) 

0 


where  for  the  flow  between  the  parallel  plates  is 


(2.9) 


y  and  r  are  measured  from  the  centerline,  ro  is  the  radius  of  the  tube  and  h  is  the  distance 
from  the  centerline  to  the  plate.  [Ref  1] 

The  energy  equations  for  the  cylindrical  and  rectangular  coordinates  are 


di  di 

ox  or 


{kr^] 

1  a 

a 

(M 

r  dr 

1 

r^d(^ 

1 

i  ar 

fV 

K  SxJ_ 

=  0 


(2.10) 


and 


di  di  di 

up - hyp - Vwp - 

dx  dy  dz 


d  ( 

dx 


dx 


V 


dy 


r 


+  - 


dz 


k^ 

V 


0  (2.11) 
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where  is  the  enthalpy.  These  equations  are  applied  to  the  problem  including  the 
boundary  conditions  and  the  temperature  distribution  is  calculated  analytically  or 
numerically.  [Ref.  2] 

2.  The  Thermally  Fully  Developed  Flow 

In  thermally  fully  developed  case,  the  relative  shape  of  the  profile  no  longer 
changes  although  the  temperature  profile  changes  with  x.  The  criterion  for  a  fully 
developed  temperature  profile  for  a  circular  tube  is 

=  const  (2.12) 

/r=Hi 

This  condition  is  reached  in  a  duct  flow  whether  the  case  is  uniform  heat  flux  (q 
through  the  wall  is  constant)  or  uniform  wall  temperature  (The  wall  temperature  is 
constant).  These  two  cases  have  a  lot  of  applications  in  engineering.  For  example, 
constant  wall  temperature  in  boiling  or  condensation,  or  constant  heat  flux  using  the 
electrical  heater.  In  the  thermally  fiilly  developed  flow,  the  convection  coefiBcient  “h”  is 
constant,  independent  of  x.  [Ref.  1] 


SrU-T. 


B.  BACKGROUND  STUDIES 


The  heat  transfer  solution  for  laminar  fully  developed  parabolic  velocity  profile 
flow  inside  a  circular  tube  in  the  thermal  region  and  subject  to  uniform  tube  wall 
temperature  was  treated  for  the  first  time  by  Graetz  in  1883  and  is  known  in  the  heat 
transfer  literature  as  Graetz  problem.  In  his  problem  Graetz  neglected  the  axial 
conduction  and  solved  the  problem.  [Ref.  3] 

The  problem  statement  is: 


1  L  AdO.  ^  1  d9, 

2  *  dx,  dxi  r*  dr. 


0,=O 

69  ^ 
dr  ~ 


at  r,  =  1 


at  r.  =  0 


(2.13) 


6,  =  I  at  X*  =  1 
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where 


e,=- . 

Tin  -To 


Re^Pr 


and  To  is  the  wall  temperature,  Tin  is  the  initial  temperature,  r  is  the  radial  distance,  x  is 
the  axial  distance,  D  is  the  diameter.  Rep  is  the  Reynolds  number  and  Pr  is  the  Prandtl 
nmnber.  [Ref.  4] 

The  velocity  profile  is  fully  developed  much  before  the  temperature  profile  when 
Ihe  Prandtl  number  is  hi^  relative  to  one.  For  such  situations,  Graetz  solution  is 
reasonably  well  justified.  But  the  Prandtl  number  is  very  small  when  the  fluid  is  liquid 
metal  and  therefore  the  axial  conduction  effects  may  not  be  neglected.  For  such  cases, 
since  the  axial  conduction  effects  are  not  included  to  the  result,  the  solution  may  give 
inaccurate  results. 

Michelsen  and  Villadsen  investigated  the  Graetz  problem  with  axial  heat 
conduction  for  circular  tube  constant  wall  temperature  case  by  using  a  numerical 
procedure.  They  used  a  method  that  is  the  combination  of  orthogonal  collocation  and 
matrix  diagonalization.  They  didn’t  include  the  wall  conduction  effects  in  the  problem. 

The  problem  may  be  defined  as: 


(1-xO 


2.66  1  d 


dy  X  dx\  dx  J  Pe  dy 


1  6^6 


(2.14) 


where 


2i?{v,)pc 


r 

x  =  -,  y. 


0  = 

T,-To 


T^-To 


The  boundary  conditions  are: 

V  — >■  -00  6  —  1 


x  =  0 


x=l 


dx 
6  =  0 


for  y  <  0 
for  y  >  0 
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where  Vz  is  the  average  fluid  velocity,  R  is  the  tube  radius,  z  is  the  axial  distance,  r  is  the 
radial  distance,  Tb  is  the  fluid  temperature  at  z  — >  -oo,  t  is  the  bulk  mean  temperature 
and  To  is  the  wall  temperature  at  z>0.  They  plotted  the  results,  Nusselt  number  vs.  y, 

heat  flux  vs.  y  and  ^  vs.  y  where  the  axial  heat  conduction  effects  are  easily  seen.  But 
they  didn’t  include  the  axial  heat  transfer  inside  the  wall,  which  increases  the  fluid 
temperature  very  much  upto  the  entrance  as  the  wall  conductivity  gets  larger.  [Ref.  5] 

X.  Yin  and  H.  H.  Bau  included  the  axial  conduction  effects  of  the  wall  to  the 
internal  flow  through  circular  tube,  by  using  two  parameters,  duct’s  outer/inner  radius 
ratio  and  fluid/wall  thermal  conductivity  ratio.  They  plotted  the  graphs  of  ^^temperature 
vs.  radius  of  the  circular  tube”  and  “Nusselt  number  vs.  Peclet  number”,  including  the 
axial  fluid  and  wall  conduction  effects.  [Ref.  6] 

C.  NUMERICAL  METHOD 

The  governing  eqiiations  are  steady  state  elliptic  partial  differential  equations.  The 
temperature  T(x,y)  throughout  the  domain  must  satisfy  the  governing  equation  and  the 
boundary  conditions  along  the  entire  boundary.  The  finite  difference  method  was  used  in 
this  numerical  study  to  express  the  governing  equations  numerically.  The  central,  forward 
and  backward  difference  expressions  used  in  this  numerical  study  for  the  grid  points  in 
figure(2)  are  as  follows  : 


Figure  (2).  The  grid  points  in  finite  difference  methods 
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The  central-difference  expressions  for  the  first  and  second  derivatives  with  error 
order  of  h^  are: 


dx  2h^ 


dy  2hy 


^0-1, j) 

dx^  (h^f 


d^t  _ 

(hyf 


(2.15a) 


(2.15b) 


(2.16a) 


(2.16b) 


2  • 

The  forward-difference  expression  with  error  order  of  h  is: 


dt  _  -t(U+2)+%J+l)  -%,j) 
dy  2hy 


2  ^ 

The  backward-difference  expression  with  error  order  of  h  is: 


hu+2)+%j+i)  [Ref.  7]  (2.18) 

^  2hy 

Substituting  these  finite-difference  expressions  into  the  governing  equations,  the  problem 
is  reduced  to  a  set  of  linear  algebraic  equations,  depending  on  the  number  of  grid  points 
used.  Usually  large  number  of  grid  points  is  desirable,  but  the  number  of  equations  to  be 
solved  becomes  too  large  in  this  case.  Gauss  Seidel  iteration  method  can  be  used  to  obtain 
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a  solution.  All  the  temperatures  of  the  grid  points  are  then  make  equal  to  the  initial 
temperature  at  the  beginning  of  the  numerical  solution.  The  temperatures  at  the  grid 
points  are  calculated  by  using  the  governing  equations  and  the  boundary  conditions.  After 
the  temperature  of  a  grid  point  is  calculated  using  the  finite  difference  equations,  this 
updated  temperature  is  used  for  the  next  grid  and  the  temperature  values  converge  to  the 
exact  temperature  a  little  more  in  every  iteration.  A  technique  called  over  relaxation  is 
used  to  speed  the  convergence  of  the  Gauss  Seidel  method  when  applied  to  the  elliptic 
partial  differential  equations.  This  technique  is: 

hiJ.new)  =  +  (1 "  ^)hj.old)  (2-19) 

w  is  known  as  the  relaxation  parameter.  For  overrelaxation  w  is  between  1  and  2;  for 
underrelaxation  between  0  and  1.  For  this  method,  overrelaxation  must  be  used  to  speed 
up  the  convergence.  [Ref.  7] 

Forsythe  and  Wasow  show  that  for  a  45  by  45  mesh  grid  the  optimum  value  of  w 
is  around  1.870.  They  also  point  out  that  when  this  relaxation  parameter  is  used,  the 
convergence  is  approximately  30  times  faster  than  the  usual  process  (w=l).  [Ref  8] 
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9  =  1^ 


to-te 


+  r  +  u  +  x/rn 
r  =  — ,  u  =  — ,  X  =- 

To  V  RePr 


u*=2(l-r^) 


The  governing  equation  becomes 


+2)69  d^9  1  d9  1  d'^t 

^  dr^  dr*  Pe^  dx* 
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(3.5) 


2.  Flow  Inside  the  Circular  Tube,  Constant  Heat  Flux  Case 

The  governing  equation  for  the  circular  duct  flow  is  (3.4) 

d^t  dt  ^  d^t  _u  dt 
dr^  r  dr  dx^  a  dx 


Now,  to  non-dimensionalize  the  problem,  assume  that 


9  = 


te-t 
ql  d/k  ’ 


+  2xld 

X  = - , 

Re^Pr 


u*  =2(\-r*^) 


The  governing  equation  becomes 


d^9  1  d9  1  d'^9 
dr*^  r*  dr*  Pe^  dx*^ 


(3.6) 


3.  Flow  Between  the  Parallel  Plates,  Constant  Wall  Temperature 


I 


Case 


The  viscous  energy  equation  in  cartesian  coordinate  system  for  the  flow  between 
the  parallel  plates  is 


di  di  di 

up - wp - \-wp - 

dx  dy  dz 


dx\  dx  j 


+  — 
dy 


dt'^ 

K - 


+  - 


dz 


(i  ^0 

k — 
dz) 


0  (3.7) 
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where  di  =  cdt  +—dP  and  is  the  radial  velocity  [Ref.l].  Then  the  equation  becomes 
P 


for  the  constant  k 


dt  dt  ^  i\ 

UpC - \rvpc - Vwpc - K\ 

dx  dy  dz 


dx 


dt 
dx  ) 


dy 


dt 

dt 


+  - 


dz 


5/ 

\dzj 


=  0  (3.8) 


The  temperature  distribution  on  z  direction  is  constant,  ^  ~  ^  laminar  and 


hydrodynamically  fully  developed,  v  =  0 ,  w  =  0 .  Equation  then  becomes 


‘  dh  dh 

=  w 

(£E? 

dx^  dy^ 

[k  j 

dx 


(3.9) 


a  =  klpc,ihe  thermal  diffusivity  of  the  fluid.  So  the  governing  equation  for  the  parallel 
plates  flow  is 

dh  dh  u  dt 


dx^  dy^  a  dx 


(3.10) 


Now,  to  non-dimensionalize  the  problem,  assume  that 


t,-t/  ^  h’  V’  RePr 


‘0  ‘e 

The  governing  equation  becomes 


1  d^0  ,  d^e  u*  de 

•  +  4- 


dx*^  dy*^  2 


(3.11) 


13 


4.  Flow  Between  the  Parallel  Plates,  Constant  Heat  Flux  Case 

The  governing  equation  for  the  parallel  plates  flow  is  (3.10) 


dh  u  dt 


dx^  dy^  a  dx 


Now,  to  non-dimensionalize  the  problem,  assume  that 


y 


,  y  =— ,  u  = — ,  X  = - 

h  V’  RePr 


+  _  2x1  Df,  +  3 


The  governing  equation  becomes 


1  d^e  \  d^e  dd 
■  +  4- 


2 


(3.12) 


5.  Effects  of  Wall  Conduction,  Flow  Inside  the  Circular  Tube, 
Constant  Wall  Temperature 

The  governing  equation  for  the  circular  duct  flow  is  (3.4) 

d^t  ^  \  dt  _u  dt 
dr^  r  dr  dx^  a  dx 

Now,  to  non-dimensionalize  the  problem,  assume  that 


i3  ^0  ^  +  r  ^  u 

0  =  — - ,  r  =  —  ,  M  =  — , 

rn  F’ 


^0  ~^e 


The  governing  equation  becomes 


d^O  1  de  1  d^0 


-+■ 


Sx:'*'  dr^  dr^  Pe^  dx* 


(3.13) 
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B.  BOUNDARY  CONDITIONS 


1.  Flow  Inside  the  Circular  Tube,  Constant  Wall  Temperature 
Case 


Figure  3.  The  boundary  conditions  of  circular  tube,  constant  wall  temperature  flow 

The  governing  equation  is  (3.5) 

/■  Add  d^e  1  de  i  d'^e 

^  ^  dr^  Pe^  dx^ 

where  the  parameters  of  the  equation  are  defined  above.  The  boundary  conditions  for  this 
case  is: 


at  r^=0, 

o 

II 

(Symmetry  boundary  condition) 

(3.14a) 

at  r'^=l, 

■0  =  0 

ISr* 

for 

for 

^  H 

+  4- 

A  IV 

o  o 

(3.14b) 

as  x'^  -00 , 

0^1 

(3.14c) 

as  x'^  ->+00 , 

0^0 

(3.14d) 
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2.  Flow  Inside  the  Circular  Tube,  Constant  Heat  Flux  Case 


Figure  4.  The  boundary  conditions  of  circular  tube,  constant  heat  flux  flow 


The  governing  equation  is  (3.6) 

L  +2]d0  8^0  1  80  1  8^0 

\l-r  I — —  = - -  + - -h - 

^  ^  8r^  Pe^ 

where  the  parameters  of  the  equation  are  defined  above.  The  boundaiy  conditions  for  this 
case  is: 


at  r^=0. 


de 

dr* 


=  0 


(Symmetry  boundary  condition) 


atr'^=l. 


dd  1 

dr*~  2 

d6 


dr* 


•  =  0 


for  >  0 


asx'"->-oo,  0^0 


as  X*  — y  +00 ,  0  — ^ 


4  2 

7  r*  r* 

■  +  ■ 


48  8  2 


(3.15a) 


(3.15b) 


(3.15c) 

(3.15d) 
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3.  Flow  Between  the  Parallel  Plates,  Constant  WaU  Temperature 

Case 


y  Insulated  wall 

^  ^  ^  yf  1  Twin  =  To  —  constant 

- ►  y*f 

^->1  ^  /  u(yl  1  : 

/  \  -  -  - .  _ 

^  ^  ^  ^  ^ '  Constant  wall  temperature 

Insulated  wall 

Figure  5.  The  boundary  conditions  of  parallel  plates,  constant  wall  temperature  flow 


The  governing  equation  is  (3.1 1) 


d'^e 


+  4 


d^0 


dx'^  dy^ 


u*  dO 
2  dx* 


where  the  parameters  of  the  equation  are  defined  above.  The  boundary  conditions  for  this 
case  is: 


at  y^=0, 

(Symmetry  boundary  condition) 

(3.16a) 

1^  =  0 

for 

for 

aty^=l, 

II 

o 

+  + 

A  IV 
o  o 

(3.16b) 

as  X* 

e^\ 

(3.16c) 

as  X*  +00 , 

e^o 

(3.16d) 
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4.  Flow  Between  the  Parallel  Plates,  Constant  Heat  Flux  Case 


Insulated  wall 

U 

_J 

u 

u 

y  y"!  q=q.-=»Baaffl 

y  Constant  heat  0ux 

Insulated  wall 

— 

n 

* 

r 

Figure  6.  The  boundary  conditions  of  parallel  plates,  constant  heat  flux  flow 


The  governing  equation  is  (3.12) 

(  1  de  1  d^e 

where  the  parameters  of  the  equation  are  defined  above.  The  boundary  conditions  for  this 
case  is; 


at  y^=0, 

at  y"'=l, 


— —  =  0  (Symmetry  boundary  condition) 

dr^ 

'  86  _  1 

dr*  4  for  X*  ^0 

—  =  0  X*  <0 

.8r* 


as  X*  ->-<10,  /  te ,  6-^0 


as  X*  ->  +00 , 


39  3 

1 

+ 

.,,.■11 

1120  16 

6  " 

(3.17a) 


(3.17b) 


(3.17c) 

(3.17d) 


18 


5.  Effects  of  Wall  Conduction,  Flow  Inside  the  Circular  Tube, 


Constant  Wall  Temperature 


Figure  7.  The  boundary  conditions  of  circular  tube,  constant  wall  temperature  including  the  circitiar  wall 
The  governing  equation  is  (3.5) 

a2. 


^  ''^+  a..+ 


80  1  8^0 


'8x^  8r^  Pe^  dx*" 

The  boundary  conditions  for  this  case  is: 


at  r^=0, 


atrM, 


atr^=l+<^,  -I 


dr* 

(Synmietiy  boundary  condition) 

(3.18a) 

^wall  ~ 

'  ^ fluid 

'  ^^wall 

_  fluid 

(3.18b) 

[  dr^ 

dr^ 

'80  _ 
<  8r^ 

0  for 

o 

VI 

+ 

(3.18c) 

0  =  0 

for 

x'">0 

where  <5r  is  — ,  t  is  the  thickness  of  the  wall,  and  r6  is  the  inner  radius  of  the  tube. 


as  X  — ^  — 00 ,  t  — ^  tgy  0  — ^  1 

as  ^  +00 ,  t  — >  Tq  ,  0  — >  0 


(3.1 8d) 
(3.18e) 
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The  Mte  difference  approximations  for  the  first  and  second  derivatives  with  error  of 

order  are  (figure  2.) 


dO 

^(.>1,;)  ^O-U) 

(4.1a) 

dx^  {ij) 

,  '^K 

dO 

+oihl) 

(4.1b) 

dr^  (ij 

)  2A, 

JjL 

(4.1c) 

A 
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2  ,2  ^\^y ) 

Sr*"  (,n  hi 


(4.1d) 


When  these  approximations  are  plugged  into  the  governing  equation  (3.5),  the  equation 
becomes 

H  Jl  li  a,  )  (4  2-1 


After  simplifying  the  expression  (4.2),  6,^ ..  becomes 


M  I  M ,  ,  M ,  ,9  r  ,  '  ,  /?  h  1 

19/4.1  n  +  "r: - tt-  +t9i_i  a  _  ..  +£'/. /..ii  tt  +  tt -  +6/1  y.n  tt-t^ — 


“"'■■'1  Pe%h!  *  2K  2h.  r  "'-•■^>1  Pe^Jit  2h,  2h  f  ‘'-^^'>1  2h,y(,,  j^)  2h^y^,_..  1  (4-3) 


Pelh]  2K  2K ) 
2  2^ 
Pelhl^  hi 


2.  The  Boundary  Conditions 


The  boimdary  conditions  are: 


a)  At  the  Center 


Ax  +  A  oa 

At  r  =0,  —  =  0 

dr^ 


(Symmetry  boundary  condition) 


1  d9 

When  this  boundary  condition  is  applied  to  the  governing  equation,  the  element - 

r*  dr* 


becomes  — .  After  applying  the  L'Hospital's  rule,  the  governing  equation  for  the  given 
boundary  condition  becomes 
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(4.4) 


d9  1  d^9 

=  2 - r  +  - 


dx*  Pe^  dx*^ 


Since  the  temperature  distribution  at  the  center  on  the  radial  direction  is  symmetric, 

<^(,•,7+1)  =  ■ 

After  applying  the  finite  difference  method,  the  at  the  boundary  becomes 


9, 


('+1,7) 


^(iJ)  - 


V  j 


1  1 

-  +  - 


\Pe\h\  '2‘hj^  j 


^(',7+1) 


^4^ 


2  '-2  ^ 
y 


\^yj 


(4.5) 


b)  At  the  Wall 


atr^=L 


^(',7)=0 


39 


=  0 


for  X*  >0 

jc^<0 


(4.6) 


FOR  xVo 


x^<0,  r‘"=l. 


39 

3r^ 


=  0 


The  governing  equation  for  the  given  boimdary  condition  becomes 

3^9  1  3^9 


Pel  dx^^ 


=  0 


(4.7) 


dO 

Since  the  wall  is  insulated,  =  0 ,  y+i)  =  0(ij-x)  at  the  boundary. 


After  applying  the  finite  difference  approximations,  the  9,^ ..  at  the  boundary  becomes 
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where  U  is  the  inlet  temperature  (as  x ->-<»)  and  to  is  the  constant  wall 
temperature.(Fluid  exit  temperature  as  x  ->  +oo  ) 

Then  the  bulk  mean  temperature  becomes 

e  ^  ^  (4.11b) 

0  is  found  out  from  the  numerical  calculation.  After  that  0^  is  calculated  by  solving  the 
equation  (4.1 1)  numerically. 


4.  The  Nusselt  Number 


The  Nusselt  number  is  (2.3) 


k 


dQ\ 

When  the  heat  transfer  equations  ql=h(T^-T^)  and  q^=-k— 


dimensionless  parameter  r^  are  plugged  in  (2.3),  the  equation  becomes 

2  de 


,  and  the 


r=R 


Nu  =  — 


6„  dr^ 


(4.12) 


r+=l 


The  finite  difference  approximation  for  the  equation  (4.12)  using  the  backward 
difference  expression  with  error  of  order  h^  becomes 


NUf..  = 
0) 


m(0 


IK 


(4.13) 


where  “j”  is  the  value  which  corresponds  to  r^— 1. 
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5.  Heat  Absorption  Upto  the  Entrance 

Assume  a  dififerential  volume,  dV  =  Iwdr  between  x  =  -oo  and  x  =  0 .  The 
energy  stored  in  this  control  volume  is 

dE  =  Cpixtdm  (4.14) 

where  At  is  the  temperature  difference  between  entrance  temperature  and  the  steady  state 
temperature  after  the  heating  starts,  At  =  t^^^  -  4  and  dm  is  the  mass  of  the  control 

volume,  dm  =  pdV  .The  energy  absorbed  in  the  control  volume  becomes 


0  ’^0 


Eabsorbed  =  ^mp  \  j(t-te)rdrdx 


(4.15) 


—CO  0 


0 

When  J  \tjdrdx  is  subtracted  from  both  sides  and  (te-to)  is  divided  to  both 

-00  0 

sides,  the  equation  (4. 1 5)  becomes 


F 

^absorbed 


0 


-00  0 


=  2;r  J  J(^  -  V)rdrdx 


(4.16) 


t  —  t  XT  T 

where  9  = - As  the  non-dimensional  parameters  of  x*  -  ^  °  and  r'^  =  —  and 

te-to  RePr 

their  first  derivatives  are  substituted  into  the  equation  (4.16),  the  equation  becomes 


^ ‘Absorbed  ^  _2^pg  J  J(^ 


(to-te)pCpro 


(4.17) 


-00  0 
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B.  FLOW  INSIDE  THE  CIRCULAR  TUBE,  CONSTANT  HEAT 
FLUX  CASE 

1.  The  Governing  Equation 

The  governing  equation  is  (3.6) 

f  1  30  1  d^6> 

1'“'' 

where 

ff-JlZL  r^=L 

g'od/k^  ro’  Re.Pr 

The  finite  difference  approximations  for  the  &st  and  second  derivatives  with  the  error  of 
order  h^  are  (4,la,b,c,d)  (figure  2.) 

dr*  (ij)  2hy 


d^0  I  ^(7^2) 

dx*^  iu) 


d^6  _  ^(»vi-i) ^0 j-i)  I  c(/i^) 

dr*^  (,j)  \ 
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When  these  approximations  are  plugged  into  the  governing  equation  (3.6),  the  equation 


becomes 


.  ^  I  %-\J)  ~  I .  (  J-l)  -  i)  I  .  1  2 


(4.18) 


After  simplifying  the  expression  (4.18),  0,^ ..  becomes 


a  ^  ,  yp-j) _ L  4./0  1  yv.j)_^  1  1  1  ^  I  1 

Vw)  p^2i,2*  Oh  'iu  "^Vu)  +^(U+1)  72+::; -  +^(,.,-n  77-r^ 


PeM  2h  2K 


1  Pm,  1 


1  1 


24  (4-19) 

2  2  ■  - - 

Pelhl'^  hi 


2,  The  Boundary  Conditions 

The  boundary  conditions  are: 


a)  At  the  Center 


At  r^=0. 


(Symmetry  boundary  condition) 


When  this  boundary  condition  is  applied  to  the  governing  equation,  the  element  — 

r*  dr* 

becomes  — .  After  applying  the  L'Hospital's  rule,  the  governing  equation  for  the  given 
boundary  condition  becomes 


-^  =  2—  1  d^0 

dx*  dr*"'''  Pe^  dx*" 


(4.20) 
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Since  the  temperature  distribution  at  the  center  on  the  radial  direction  is  symmetric. 


After  applying  the  finite  difference  method,  the  at  the  boundary  becomes 


e, 


{.MJ) 


1  1 


yPejyhx  2hx  ^ 


1  1 


2hx 


>1 

f  4  1 

J 

K^y  ) 

2  2^ 

Pelhl  ^  hi 


(4.21) 


b)  At  the  Wall 


atr^=l. 


FOR  >  0, 


de 

dr* 


2 


k^  =  ql,  G{x*,r*)  =  - 
dr  , 


for 

X*  >0 

for 

<0 

{x,r) 

r*-- 

l/k  ’ 

fo 

(4.22) 


so  the  equation  becomes  . 

dr*  2 


The  finite  difference  approximation  for  this  equation  is 


2h^  2 


When  the  boundary  condition  is  applied  to  the  governing  equation  (3.6),  the  equation 
becomes 


d^e  1  1  ^ 

+  — :r — r  =  0 


dr*  2r*  Pe^  dx* 


(4.23) 
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As  the  finite  difference  approximations  are  substituted  in  the  equation  above,  the  equation 
becomes 


9. 


- ' 


'_J_' 


^  1  ^ 


^  ^ 


+  9, 


J__  ] 
hy  I 


-  +  ■ 


Pelhl  '  hi 


FOR  x^<0 


x^<0,  r''=l, 


39 

dr^ 


=  0 


The  governing  equation  for  the  given  boundary  condition  becomes 


3^9  1  3^9  ^ 

+  — ^ - :^  =  0 


,2  p„2  2 

3r^  dx* 


dO 

Since  the  wall  is  insulated,  =  0 ,  =  9^^  J_^■^  at  the  boundary. 

After  applying  the  finite  difference  approximations  to  (4.25),  the 
becomes 


KP^lhlj 


+  9, 


(i-i.y) 


^  1  ^ 


^2^ 


J 


2  2 

+  ■ 


Pelhl  hi 


c)  As  X  ^  -  00 


as  X*  — >  —00 ,  9^^  ^  0 


(4.24) 


(4.25) 


at  the  boundary 


(4.26) 


(4.27) 
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d)  As  X  “^  +  CO 


fls  x"*"  +00  j  ^  2x  ^“ 


7 

-  +  - 


48  8 


(4.28) 


3.  The  Bulk  Mean  Temperature 

The  mean  temperature  at  a  position  “i”  is  (2.8) 


2 

f>n=r^  l^trdr 
^^0  0 


fo 

where  V  is  the  average  velocity  and  ro  is  the  radius  of  the  tube.  When  jtjilm'dr  is 

0 

subtracted  from  both  sides  and  dimensionless  parameter  r'*’  is  substituted  inside  the 
equation,  the  equation  becomes 

-0=  !(f-  Qu^^r^dr^  (4.29) 

0 

Let’s  define  the  bulk  mean  temperature, 

(4-30a) 


where  tg  is  the  inlet  temperature  (as  x  ^  -oo ). 
Then  the  bulk  mean  temperature  becomes 


0 


(4.30b) 
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0  is  found  out  from  the  numerical  calculation.  After  that  0^^  is  calculated  by  solving  the 
equation  above  numerically. 


4.  The  Nusselt  Number 

hei 

The  Nusselt  number:  Aw  =  — . 

k 

The  heat  transfer  from  the  wall:  a"  = /zfr  ,  ^  -t  ) 

The  temperature  at  the  position  of  (x^/=l):  =t,-eq„  d/k 

The  average  temperature  at  (x"^):  9^q\  djk 

When  these  equations  are  substituted  into  the  Nusselt  number  equation,  Nusselt 
number  becomes: 


(4.31) 

5.  Heat  Absorption  Upto  the  Entrance 

Assume  a  differential  volume,  dV  =  27jrdr  between  x  =  — oo  and  x  =  0 .  The 
energy  stored  in  this  control  volume  is  (4. 1 4) 

dE  =  CpAtdm 

where  At  is  the  temperature  difference  between  entrance  temperature  and  the  steady  state 
temperature  after  the  heating  starts.  At  =  -4  and  dm  is  the  mass  of  the  control 

volume,  dm  =  pdV  .The  energy  absorbed  in  the  control  volume  becomes 
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Eabsorbed  =  '^^f^P  J  \{t  -  te)rdrdx . 


(4.32) 


When  both  sides  are  divided  by  djk ,  the  equation  (4.32)  becomes 


TT  0 

f  absorbed^  _  J  ^fffdrdx 

ioP^pfol^  -“0 


(4.33) 


x!  f* 

where  0=  As  the  non-dimensional  parameters  of  =  °  and  =  ~and 

qodik 

their  first  derivatives  are  substituted  into  the  equation  (4.33),  the  equation  becomes 


^absorbed_  ^  J  dv'^ dx^ 

(loIXp^tl^  -“0 


(4.34) 


C.  FLOW  BETWEEN  THE  PARALLEL  PLATES,  CONSTANT 
WALL  TEMPERATURE  CASE 
1.  The  Governing  Equation 

The  governing  equation  is  defined  in  equation  (3.11)  as 

1 


where 


^0 


h  RePr 


The  finite  difference  approximations  for  the  first  and  second  derivatives  with  error  of 
order  h^  are  (4.1a,b,c,d)  (figure  2.) 
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d9 

dx  (ij) 


+  o(hl) 


de 

dr*  (ij) 


d^e 

dx  (/,;•) 


^(<+1,7) 


h 


2 

X 


+  o{hl) 


d^0 

dr*"  ii,n 


+  o{hl) 


When  these  approximations  and  are  plugged  into  the  governing  equation  (3.11),  6^^ 
becomes 


2.  The  Boundary  Conditions 

The  boxmdary  conditions  are: 

a)  At  the  Center 

dG 

At  y^=0,  j=l,  — -  =  0  (Symmetry  boundaiy  condition) 

dy 
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When  this  boundary  condition  is  applied  to  the  governing  equation,  the  governing 
equation  becomes 


1  ^  d'^e  _3  39 

Pe\  dx*'  4  dx^ 


(4.36) 


Since  the  temperature  distribution  at  the  center  on  the  radial  direction  is  symmetric, 


at  the  center. 

After  applying  the  finite  difference  method  to  (4.36),  the  <9-.  .  at  the  boundary  becomes 


(1+1,1) 


yPelhl  %\j 


^  ^(i-i.i; 


1  3 

-  +  - 


\P^D^x  J 


+  0, 


^8^ 


('■.2) 


2 

Pelhl  hi 


(4.37) 


b)  At  the  Wall 


atyM,  ( 


d0 


=  0 


for  X*  >0 


for  x""  <  0 


(4.38) 


FOR  xVQ 


X  <0,  y=l,  =  0 


The  governing  equation  for  the  given  boundary  condition  becomes 

1  d^9  ,  d^e  ^ 

-  +  4 — r  =  0 


Pe^^  dx*  dy^ 


(4.39) 
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do 

Since  the  wall  is  insulated,  =  0 ,  at  the  boundary. 


After  applying  the  finite  difference  approximations,  the  at  the  boundary  becomes 


0. 


(MJ) 


^  1  ^ 


7) 


^8^ 


2  8 
Pelhl  ■"  hi 


(4.40) 


e)  As  X  -  00 

as  x""  ->-00,  ^1  (4.41) 

f)  As  X  ■>  +  00 

as;c^->+oo,  (4.42) 


3.  The  Bulk  Mean  Temperature 

The  mean  temperature  is  defined  in  equation  (2.9)  as 

When  the  dimensionless  parameters  u"^,  y"^  are  substituted  in  (2.9),  the  tm  becomes 


tm  =  \u^tdy* 

0 

where  V  is  the  average  velocity  and  h  is  the  distance  fi'om  the  centerline  to  the  plate. 
1 

When  tg  =  ItgU^dy'^  is  subtracted  from  both  sides,  the  equation  becomes 

O' 
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0 


Let’s  define  the  bulk  mean  temperature, 


a  (^o 


(4.43a) 


where  te  is  the  inlet  temperature  (as  x-^-oo)  and  to  is  the  constant  wall 
temperature.(Fluid  exit  temperature  as  x  — +oo  ) 

Then  the  bulk  mean  temperature  6^  becomes 


(4-43b) 

^  0 

6  is  taken  from  the  numerical  calculation.  After  that  6^  is  calculated  by  solving  the 
equation  (4.43)  numerically. 


4.  The  Nusselt  Number 

The  Nusselt  number  is  defined  as  (2.3) 

Nu  =  — 
k 


•I 

When  the  heat  transfer  equations  q"^  =  h(Tg  -T^),  gl  =  k-—\ 

dy 


parameter  r"^  are  plugged  in  (2.3),  the  equation  becomes 


,  and  the  dimensionless 


y^h 


Nu  =  - 


4  dd 


(4.44) 
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The  jSnite  difference  approximation  for  the  equation  (4.44)  using  the  backward  difference 
expression  with  error  of  order  h  becomes 


NUi  = 


9, 


mi 


2K 


(4.45) 


5.  Heat  Absorption  Upto  the  Entrance 

Assume  a  differential  volume,  dV  =  bdxdy  between  x  =  -oo  and  x  =  0 .  The 
energy  stored  in  this  control  volume  is  (4.14) 

dE  =  CpAtdm 

where  Aris  the  temperature  difference  between  entrance  temperature  and  the  steady  state 
temperature  after  the  heating  starts,  At  =  and  dm  is  the  mass  of  the  control 

volume,  dm  =  pdV  .The  energy  absorbed  in  the  control  volume  becomes 


■^absorbed 


0  +h 

=  bpCp  J  \(t-Qdydx. 


(4.46) 


0  h 

When  IbpCp  j  \tgdydx  is  subtracted  from  both  sides  and  {tg-to)  is  divided  to  both 

-00  0 


sides,  the  equation  (4.46)  becomes 


F 

■^absorbed 

bpcpitg-t;) 


=  2  I  ]{0-\)dydx 
—00  0 


(4.47) 
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where  0  =  -^—^  As  the  non-dimensional  parameters  of 

(4-U 


r*  =  —  anA  their  first  derivatives  are  substituted  into  the  equation 
becomes 


and 

RePr 

(4.47),  the  equation 


0  1 


ahsorbed - =  _4Pe  J  J(^  -  \)dy*dx^ 


(4.48) 


IfbpcAt^  -t,) 


-00  0 


D.  FLOW  BETWEEN  THE  PARALLEL  PLATES,  CONSTANT 
HEAT  FLUX  CASE 


1.  The  Governing  Equation 

The  governing  equation  is  (3.12) 

1  d^e  .  d^9 

- 2 - 7  +  ^ - 2 

dx^  dy^ 


u*  dO 
2  dx* 


where 


te-t 


.+  2x1  Dt, 

^  ^ 
RePr 


Dh=4h 


The  fiinite  difference  approximations  for  the  first  and  second  derivatives  are  (4.1a,b,c,d) 


at*  (.■.;) 

as  +<,(*;) 

(ij) 
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dx^ 


(i,j) 


hi 


+  oihl) 


_  ^(, ■,;+!)  +  ^(lJ-1) 


dr* 


+  o(hl) 


iU) 


When  these  approximations  and  u"*^  are  plugged  into  the  governing  equation,  the  bulk 
temperature  at  the  position  (i  j)  is  written  as 


■'(.M.j) 


^(V) 


1  ,  Hj,  3  ] 

{Pelhl  8/2,  Zh,] 

[Pelhl  ZK  8aJ 

(4.49) 


2  ^ 
Pelh^  hi 


2.  The  Boundary  Conditions 

The  boundary  conditions  are: 


a)  At  the  Center 


At  y^-0,  j=l,  — 7-^  (Symmetry  boundary  condition) 

dy 

When  this  boundary  condition  is  applied  to  the  governing  equation  (3.12),  the  governing 
equation  becomes 


1  ^  ^  d^e  _3  dO 

Pel^dx*"^  dy*"~Adx* 


(4.50) 


Since  the  temperature  distribution  at  the  center  on  the  radial  direction  is  symmetric, 
at  the  center. 
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After  applying  the  finite  difference  method,  the  at  the  boimdary  becomes 


d. 


(z+1.1) 


yPelhl 


1  3 

+  ■ 


yPelhl  ik 


2  8 


P^dK  K 


fsl 

1^2 

J 

J 

(4.51) 


,  +  1  d9 

aty=i,  —  = 

dy 


b)  At  the  Wall 

0 


4 


for 

<0 

for 

+ 

IV 

o 

FOR  x^<0 


x^<0,  y'^=l,  =  0 


The  governing  equation  for  the  given  boundary  condition  becomes 

1  d^0  ,  d^0  „ 

■  +  4 — ;-  =  0 


P^\  dx*  dy*^ 


(4.52) 


d0 


Since  the  wall  is  insulated,  ^ at  the  boundary. 


After  applying  the  finite  difference  approximations,  the  0^^  ^^  at  the  boundary  becomes 


^iij)  ~  ' 


yP^lhl 


r 


K 


2 

Pelhl  ^  hi 


(4.53) 
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FOR  >  0 


The  governing  equation  for  the  given  boundary  condition  becomes 

1  d^0  ,  d^e  ^ 

Peo.  dx*  dy* 


(4.54) 


dy*  4  ’  2h.  4 


4  ’  -  y  boundaiy. 


After  applying  the  finite  difference  approximations,  the  at  the  boimdary  becomes 


^(<  j  )  “ 


PelK 


2  7.2  ■'■^(>■-1.;)  d^2t2 


1  Li 

Pelhl) 


(4.55) 


PelK  K 


c)  As  X  -  00 


as  -CO ,  P(ij)  ^ 


(4.56) 


d)  As  X  ->  +  00 


as  +CO ,  6,,  =  -2x  + 


+  .  39  ,3  y 


1120  16  6 


(4.57) 


3.  The  Bulk  Mean  Temperature 


The  mean  temperature  is  (2.9) 
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When  the  dimensionless  parameters  u"^,  are  substituted  in  (2.9),  the  tm  becomes 


r,  =  'ju*Tc!y* 

0 

where  V  is  the  average  velocity  and  h  is  the  distance  from  the  centerline  to  the  plate. 
1 

When  is  subtracted  from  both  sides,  the  equation  becomes 

0 

0 

Let’s  define  the  bulk  mean  temperature, 


n  _  \  £ 


(4.58a) 


where  te  is  the  inlet  temperature  (as  x  — >  -oo ).  Then  the  bulk  mean  temperature  9„ 


becomes 


9. 


(4.58b) 


9  is  taken  fi:om  the  numerical  calculation.  After  that  9j^  is  calculated  by  solving  the 
equation  (4.58)  numerically. 


4.  The  Nusselt  Number 

The  Nusselt  number:  JVw  =  . 

k 

The  heat  transfer  from  the  wall:  q'^  =  )  • 

The  temperature  at  the  position  of  (xV^=l):  t  *  *  ■=t^-9q"^Df,lk 

(JT  —1) 
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The  average  temperature  at  (x"^:  -  6„ql  D^/k 

When  these  equations  are  substituted  into  the  Nusselt  number  equation,  Nusselt  number 
is  calculated  as 


Nu  =  - 


9„-e 


(4.59) 


5.  Heat  Absorption  Upto  the  Entrance 

Assume  a  differential  volume,  dV  =  hdxdy  between  x  =  -oo  and  x  =  0 .  The 
energy  stored  in  this  control  volume  is  (4.14) 

dE  =  CpAtdm 

where  At  is  the  temperature  difference  between  entrance  temperature  and  the  steady  state 
temperature  after  the  heating  starts,  M  =  -t^  and  dm  is  the  mass  of  the  control 

volume,  dm  =  pdV  .The  energy  absorbed  in  the  control  volume  becomes 


0 

^ absorbed  ~  ^ p  J 

-00  —h 


(4.60) 


When  both  sides  are  divided  by  the  equation  (4.60)  becomes 


absorbed  = -g  J 

qj)pcphlk  _„o 


(4.61) 


where  6  =  .  As  the  non-dimensional  parameters  of  x'*’  =  — and  =  —  and 

qoDjk  RePr  r, 

their  first  derivatives  are  substituted  into  the  equation  (4.61),  the  equation  becomes 


=  -\6Pe  I  \edy-^dx^ 

qJbhVcc 


(4.62) 
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E.  THE  EFFECTS  OF  WALL  CONDUCTION,  FLOW  INSIDE  THE 


CIRCULAR  TUBE,  CONSTANT  WALL  TEMPERATURE  CASE 
1.  The  Governing  Equation 

The  governing  eqxmtion  is  (3.5) 

^2)  de  d^d  1  dd  1  d^e 

!-'■  — r  =  — T  +  — r— r  +  — 5 - T 

I  Jdx  5^+2  ^  Qr  Pen  0^*^ 

where 

to-te  ro  RePr 


The  finite  difference  approximations  for  the  first  and  second  derivatives  with  error  of 
order  h^  are  (4.1a,b,c,d)  (figure  2.) 
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d^e 


dr^ 


(ij) 


hi 


+  o{hl) 


When  these  approximations  ate  plugged  into  the  governing  equation  (3.5),  the  equation 
becomes 


Pei 


ij) 


hi 


ynj) 


2/2. 


-(l- 


y^j) 


,2  1  I  _ 


2K 


=  0 


(4.63) 


After  simplifying  the  expression  (4.63),  the  equation  becomes 


~  ■ 


^  I  yl.n  _L 

^Pe]^l  2K  IK 


_! _ 

Pelhl  2\  Ih, 


+  9, 


vjm 


1  1 


\K  ^Kpo-dj 


+9, 


('./-O 


J _ 1_ 


(4.64) 


Pelhl  hi 


2.  The  Boundary  Conditions 

The  boundary  conditions  are: 


a)  At  the  Center 

+  dO 

At  r^=0,  —  =  0  (Symmetiy  boundary  condition) 

dr* 


1  39 

When  this  boundary  condition  is  applied  to  the  governing  equation,  the  element - 

r*  dr* 

becomes  ^ .  After  applying  the  LUospital's  rule,  the  governing  equation  for  the  given 
boundary  condition  becomes 


39  ^ 


3x*  dr*^  Pe^  dx*' 


(4.65) 


Since  the  temperature  distribution  at  the  center  on  the  radial  direction  is  symmetric. 
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After  applying  the  finite  difference  method,  the  Of^  j^  at  the  boundary  becomes 


^(/j)  “ 


^Pelhl  2Kj 


+  ^(M  j) 


"  1  ^  1  ^ 
sJPelhl  ^  2h^j 


2  _4_ 

Pelhl  ^  /zj 


b)  At  the  Wall 


at  r^=l, 


^wall  ^ fluid 

derail  _  fluid 


[  dr^  dr^ 

The  heat  fluxes  through  the  fluid  and  the  wall  are  equal  to  each  other  at  r^ 

T.  ^^solid  _  u  fluid 

dr^  dr'^ 

Using  the  2“**  order  forward  and  backward  difference  method,  the 
becomes 


V 


~  %j-2)  ^ 

^  '^^yf  J 

-^s 

If  hys=hyf,  the  bulk  temperature  at  the  fluid- wall  boundary  becomes 


r 


%j+i) 

.  _ I 


w  kf  ^ 

<9(,j+2)+-^%y-2) 


1  +  ^ 

V  J 


(4.66) 


1. 


equality 


(4.67) 


(4.68) 
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c)  On  the  Outer  Edge  of  the  Wall 


FOR  x^<0 


For  x'^<0  and  =  1  +  ,  the  outer  surface  of  the  solid  wall  is  insulated  where 


solid  _  wall  _  thickness 
tube  inner  radius 


The  governing  equation  for  the  given  boundary  condition  becomes 

d^9  1  d^e  ^ 

■  +  — ^ — r  =  0 


Pe^  dx^ 


(4.69) 


d9 

Since  the  wall  is  insulated,  — -  =  0 ,  =  9.^  at  the  boimdary. 

dr* 

After  applying  the  finite  difference  approximations,  the  9q  j-^  at  the  boundary  becomes 


^(MJ) 


^  1  ^ 


+  ^(1-1,/ 


"  1  ^ 


^2^ 


2  2 

■  +  ■ 


Pw  K 


(4.70) 


FORx^>l 


%)=0 


(4.71) 


d)  As  X  -  00 

as  a:*  ->  -00 ,  ^(ij)  f  (4-72) 

e)  As  X  “^  +  00 
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as  ->  +05 , 


^(.u)  ^ 


(4.73) 


3.  The  Bulk  Mean  Temperature 

The  mean  temperature  at  a  position  “i”  is  (2.8) 

2 

^m(x)  ~  jr  2 

^^0  0 

ro 

where  V  is  the  average  velocity  and  r©  is  the  radius  of  the  tube.  When  ^t^ulm-dr  is 

0 

subtracted  from  both  sides  and  dimensionless  parameter  r"^  is  substituted  inside  the 
equation,  the  equation  becomes 

0 


Let’s  define  the  bulk  mean  temperature, 


,  where  te  is  the  inlet  temperature 


(as  X  ->  -00  )  and  to  is  the  constant  wall  temperature.(Fluid  exit  temperature  as  x  ->  +oo ) 
Then  the  bulk  mean  temperature  becomes 


(4.74) 


0 

9  is  found  out  from  the  numerical  calculation.  After  that  9^  is  calculated  by  solving  the 
equation  (4.74)  nmnerically. 
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4.  The  Nusselt  Number 


The  Nusselt  number  is  dejSned  as  (2.3) 


Ar  ^ 

Nu  =  — . 
k 


When  the  heat  transfer  equations  qo=Wo~K)  ^nd  q"o=-k 


dimensionless  parameter  r"^  are  plugged  in  (2.3),  the  equation  becomes 


dr 


,  and  the 


r=R 


Nu  =  - 


2  dO 


dr'' 


r+=i 


(4.75) 


The  finite  difference  approximation  for  the  equation  (4.75)  using  the  backward 
difference  expression  with  error  of  order  h^  becomes 


Nun.  = - 


2h 


(4.76) 


where  “j”  is  the  value  which  corresponds  to  r^=l. 


5.  Heat  Absorption  Upto  the  Entrance 

Assume  a  differential  volume,  dV  -  27irdr  between  x  =  -oo  and  x  =  0 .  The 
energy  stored  in  this  control  volume  is  (4.14) 

dE  =  CpAtdm 

where  At  is  the  temperature  difference  between  entrance  temperature  and  the  steady  state 
temperature  after  the  heating  starts.  At  =  -4  and  dm  is  the  mass  of  the  control 

volume,  dm  =  pdV  .The  energy  absorbed  in  the  control  volume  becomes 
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Eabsorbed  =  '^^P^P  \  W-te)rdrdx.  (4.77) 

—00  0 


0  ''o 

When  iT^pCp  J  jt^rdrdx  is  subtracted  from  both  sides  and  (te-to)  is  divided  to  both 

—00  0 


sides,  the  equation  (4.77)  becomes 


P  0 

■  =2^  I  ](&  -  \)rdrdx 

ife  — «>  0 


(4.78) 


^ ^  X  f  y 

where  9  = - — .  As  the  non-dimensional  parameters  of  ^  °  and  r'^  =  —  and 

t  -t  RePr  r. 


their  first  derivatives  are  substituted  into  the  equation  (4.16),  the  equation  becomes 


r  0  1 

— ^^bsorb^  ^  _2^Pe  I  _  \)r*dr*dx^  (4.79) 

{to -te)  pc  pro  -®o 
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V.  NUMERICAL  RESULTS  AND  DISCUSSION 


A.  BULK  MEAN  TEMPERATURE  VS.  X 

The  bulk  mean  temperature  graphs  of  the  fully  developed  laminar  flow  inside  the 
circular  tube  and  parallel  plates  for  constant  wall  temperature  and  heat  flux  cases  were 
plotted  in  figures  8-11.  The  Peclet  numbers  0.5,  1,  2,  3,  5,  10,  50,  and  100  were  used  for 
the  plots  to  see  how  the  Peclet  number,  Pe=RePr,  affects  the  axial  conduction  in  the  flow. 
As  it  is  seen  from  the  figures,  when  the  Peclet  number  gets  smaller,  the  axial  conduction 
effects  of  the  fluid  becomes  very  important  and  the  temperature  of  the  fluid  upto  the 
entrance  increases.  The  axial  conduction  of  the  fluid  has  to  be  considered  for  these  cases 
and  must  not  be  neglected  for  the  Peclet  numbers  less  than  about  10.  When  the  Peclet 
number  gets  larger,  the  effects  of  the  axial  conduction  decrease.  As  it  is  seen  from  the 
figures  8-11,  for  the  Peclet  numbers  50  and  100,  the  bulk  mean  temperature  of  the  fluid  at 
the  position  x‘'=0  is  almost  equal  to  the  entrance  temperature,  which  means  that  there  is 
almost  no  heat  absorption  or  temperature  increase  upto  the  entrance.  The  axial  conduction 
of  the  fluids  can  be  neglected  for  these  high  Peclet  numbers. 

1.  Flow  Inside  the  Circular  Tube,  Constant  Wall  Temperature 

Case 

The  graphs  of  the  bulk  mean  temperature  vs.  x  for  the  circular  tube,  constant  wall 
temperature  flow  was  plotted  for  various  Peclet  numbers  in  figure  8. a  and  8.b  where  theta 
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bulk  mean  is 


"■fc-0 


and  the  nondimensional  axial  distance  x'*'  is  x* 


2x/d 
RePr  ' 


The  validity  of  the  figure  8.b  has  been  checked  by  comparison  with  Ref.5. 


CIRCULAR  TUBE  CONSTANT  WALL  TEMP.  FLOW,  THETA  BULK  MEAN 


Figure  8.a  The  axial  bulk  mean  temperature  distribution  for  circular  tube  constant  wall  temperature  flow 


CIRCULAR  TUBE  CONSTANT  WALL  TEMP.  FLOW.  THETA  BULK  MEAN 


Figure  8.b  The  enlarged  view  of  figure  8.a 
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2.  Flow  Inside  the  Circular  Tube,  Constant  Heat  Flux  Case 

The  graphs  of  the  bulk  mean  temperature  vs.  x  for  the  circular  tube,  constant  heat 
flux  flow  was  plotted  for  various  Peclet  numbers  in  figure  9.a  and  9.b  where  theta  bulk 


mean  is  0„ 


qodik 


and  the  nondimensional  axial  distance  x"^  is 


2x1  d 
Re^  Pr  ■ 


CIRCULAR  TUBE  CONSTANT  HEAT  FLUX  FLOW,  THETA  BULK  MEAN 


Figure  9.a  The  axial  bulk  mean  temperature  distribution  for  circular  tube  constant  heat  fltix  flow 


CIRCULAR  TUBE  CONSTANT  HEAT  FLUX  FLOW.  THETA  BULK  MEAN 


Figure  9.b  The  enlarged  view  of  figure  9.a 
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3.  Flow  Between  the  Parallel  Plates,  Constant  Wall  Temperature 
Case 

The  graphs  of  the  bulk  mean  temperature  vs.  for  the  parallel  plates,  constant  wall 
temperature  flow  was  plotted  for  various  Peclet  numbers  in  figure  lO.a  and  lO.b  where 

theta  bulk  mean  is  ^  and  the  nondimensional  axial  distance  x""  is  jc+  =  . 

"  {to-te)  RePr 


PARALLEL  PLATES  CONSTANT  WALL  TEMP.  PLOW,  THETA  BULK  MEAN 


X 


Figure  lO.a  The  axial  bulk  mean  temperature  distribution  for  parallel  plates  constant  wall  temperature  flow 


PARALLEL  PLATES  CONSTANT  WALL  TEMP.  FLOW.  THETA  BULK  MEAN 


X 


Figure  lO.b  The  enlarged  view  of  figure  lO.a 


4.  Flow  Between  the  Parallel  Plates,  Constant  Heat  Flux  Case 

The  graphs  of  the  bulk  mean  temperature  vs.  x  for  the  parallel  plates,  constant  heat 
flux  flow  was  plotted  for  various  Peclet  numbers  in  figure  1 1  .a  and  1 1  .b  where  theta  bulk 


mean  is  6^^ 


(I'oDjk 


and  the  nondimensional  axial  distance  x"^  is  = 


2x1  Df, 

RePr  ‘ 


PARALLEL  PLATES  CONSTANT  HEAT  FLUX  FLOW,  THETA  BULK  MEAN 


Figure  ll.a  The  axial  bulk  mean  temperature  distribution  for  parallel  plates  constant  heat  flux  flow 

PARALLEL  PLATES  CONSTANT  HEAT  FLUX  FLOW.  THETA  BULK  MEAN 


Figure  ll.b  The  enlarged  view  of  figure  1  La 
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B.  NUSSELT  NUMBER  VS.  X 


The  nusselt  number  graph  of  the  fully  developed  laminar  flow  inside  the  circular 
tube  and  parallel  plates  for  constant  wall  temperature  and  heat  flux  cases  were  plotted  in 
figures  12-15.  The  Peclet  numbers  0.5,  1,  2,  3,  5,  10,  50,  and  100  were  used  for  the  plots 
to  see  how  the  Peclet  number,  Pe=RePr,  affects  the  nusselt  number  distribution  in  the 
flow. 

1.  Flow  Inside  the  Circular  Tube,  Constant  Wall  Temperature 

Case 

The  graph  of  the  nusselt  number  vs.  x  for  the  circular  tube,  constant  wall 
temperature  flow  was  plotted  for  various  Peclet  numbers  in  figure  12  where  x"^  is 
+  2x1  d 

X  = - . 

RePr 


CIRCULAR  TUBE  CONSTANT  WALL  TEMP.  FLOW.  NUSSELT  NUMBER 


X 


Figure  12.  The  nusselt  number  distribution  for  circular  tube  constant  wall  temperature  flow 
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The  nusselt  numbers  approach  to  the  values  in  those  of  table  1  when  x  gets  larger  as  seen 


in  figure  12. 


Pe 

0.5  1  2  3  5 

10  50  100 

Nu.  no 

4.0971  4.0276  3.9224  3.8506  3.767 

3.6948  3.6586  3.6572 

Table  1.  T 

tie  nusselt  number  for  various  Peclet  numbers  for  circular  tube  constant  wall  temperature  case  as 

X  goes  to  infinity 


2.  Flow  Inside  the  Circular  Tube,  Constant  Heat  Flux  Case 


The  graph  of  the  nusselt  number  vs.  x  for  the  circular  tube,  constant  heat  flux  flow 


was  plotted  for  various  Peclet  numbers  in  figure  13  where  x  is  x  =  — — — . 

Re^  Pr 


CIRCULAR  TUBE  CONSTANT  WALL  TEMP.  FLOW,  NUSSELT  NUMBER 


Figure  13.  The  nusselt  number  distribution  for  circular  tube  constant  heat  flux  flow 


The  nusselt  number  is  independent  from  the  Peclet  number  for  this  case,  and  approaches 
to  4.364  where  x  goes  to  infinity. 
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3.  Flow  Between  the  Parallel  Plates,  Constant  Wall  Temperature 


Case 


The  graph  of  the  nusselt  number  vs.  x  for  the  parallel  plates,  constant  wall 
temperature  flow  was  plotted  for  various  Peclet  numbers  in  figure  14  where  is 


2x1  Df^ 
RePr 


,  is  4h. 


PARALLEL  PLATES  CONSTANT  WALL  TEMP.  FLOW.  NUSSELT  NUMBER 


Figure  14,  The  nusselt  number  distribution  for  parallel  plates  constant  wall  temperature  flow 


The  nusselt  numbers  approach  to  the  values  in  those  of  table  2  when  x  gets  larger  as  seen 
in  figure  12. 


Pe 

1 

2 

3 

5 

10 

50 

100 

Nu.  no 

8.0588 

8.0059 

7.9168 

7.8464 

7.7471 

7.6303 

7.5457 

7.542 

Table  2.  The  nusselt  number  for  various  Peclet  numbers  for  parallel  plates  constant  wall  temperature  case 
as  X  goes  to  infinity 
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Flow  Between  the  Parallel  Plates,  Constant  Heat  Flux  Case 


The  graph  of  the  nusselt  number  vs.  x  for  the  parallel  plates,  constant  heat  flux 


flow  was  plotted  for  various  Peclet  numbers  in  figure  15  where  x""  is  x  = 


.  .  ,  2x1  Du  ^  . 

\DhlS 


RePr 


4h. 


PARALLEL  PLATES  CONSTANT  HEAT  FLUX  FLOW,  NUSSELT  NUMBER 


Figure  15.  The  nusselt  number  distribution  for  parallel  plates  constant  heat  flux  flow 


The  nusselt  number  is  independent  from  the  Peclet  number  for  this  case,  and 
approaches  to  8.235  where  x  goes  to  infinity. 

C.  HEAT  ABSORPTION  UPTO  THE  ENTRANCE 

The  energy  change  from  the  entrance  to  x=0  where  the  heating  starts  has  been 
investigated  and  plotted  for  fully  developed  laminar  flow  inside  the  circular  tube  and  parallel 
plates  for  constant  wall  temperature  and  heat  flux  cases  in  figures  16-19.  The  Peclet  numbers 
0.5,  1,  2,  3,  5,  10,  50,  and  100  were  used  for  the  plots  to  see  how  the  Peclet  number,  Pe=RePr, 
affects  the  heat  absorption  of  the  fluid.  As  it  is  seen  from  the  graphs,  the  heat  absorption 
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logarithmically  increases  as  the  Peclet  number  decreases.  Therefore  axial  conduction  effects  or 
heat  absorption  upto  the  entrance  where  the  heating  starts  has  to  be  considered  for  low  Peclet 


numbers. 


1.  Flow  Inside  the  Circular  Tube,  Constant  Wall  Temperature 


Case 


The  graph  of  the  heat  absorption,  E*  vs.  Peclet  number  for  the  circular  tube, 
constant  wall  temperature  flow  was  plotted  in  figure  16  where  the  dimensionless 


absorbed  energy  E*  is  £*  ^ ^ 


ENERGY  ABSORPTION,  CIRCULAR  TUBE,  CONSTANT  WALL  TEMP. 


10"  10’ 
Peclet  Number 


Figure  16.  The  energy  absoption  upto  the  entrance  for  circular  tube  constant  wall  temperature  flow 


2.  Flow  Inside  the  Circular  Tube,  Constant  Heat  Flux  Case 

The  graph  of  the  heat  absorption,  E*  vs.  Peclet  number  for  the  circular  tube, 
constant  heat  flux  flow  was  plotted  in  figure  17  where  the  dimensionless  absorbed  energy 


IS 
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Figure  17.  The  energy  absoption  upto  the  entrance  for  circular  tube  constant  heat  flux  flow 


3.  Flow  Between  the  Parallel  Plates,  Constant  Wall  Temperature 

Case 

The  graph  of  the  heat  absorption,  E*  vs.  Peclet  number  for  the  parallel  plates,  constant 
wall  temperature  flow  was  plotted  in  figure  18  where  the  dimensionless  absorbed  energy  E  is 


r  =• 


^absorbed 


hrbpcp{t„  -tj 


Figure  18.  The  energy  absoption  upto  the  entrance  for  parallel  plates  constant  wall  temperature  flow 
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4.  Flow  Between  the  Parallel  Plates,  Constant  Heat  Flux  Case 

The  graph  of  the  heat  absorption,  E*  vs.  Peclet  number  for  the  parallel  plates,  constant 
heat  flux  flow  was  plotted  in  figure  19  where  the  dimensionless  absorbed  energy  E*  is 


^absorbed 

glbh^a' 


Figure  19.  The  energy  absorption  upto  the  entrance  for  parallel  plates  constant  heat  flux  flow 


D.  AXIAL  CONDUCTION  EFFECTS  INCLUDING  THE  WALL 
CONDUCTION 

The  axial  conduction  effects  for  circular  tube  constant  wall  temperature  case  has 
been  investigated  including  the  wall  conduction  effects.  The  Peclet  numbers  1,  3,  5, 
outer/inner  radius  ratios  of  1.2, 2  and  fluid/wall  thermal  conductivity  ratios  of  1  and  1/10 
were  used  to  see  how  these  properties  affect  the  temperature,  nusselt  number  and  heat 
absorption.  The  theta  buUc  mean  temperature,  nusselt  nmnber  variations  were  plotted  vs. 


64 


x""  and  the  heat  absorption  were  quantified  as  a  function  of  Peclet  number,  thermal 
conductivity  and  radius  ratios.  The  results  are  seen  on  the  following  pages. 

1.  Bulk  Mean  Temperature  Variations 


The  vs.  graphs  are  seen  in  figures  20.a  and  20.b  for  different  Peclet 

^cl,x=0 

numbers,  thermal  conductivity  and  wall  thickness.  Since  the  bulk  mean  temperature 
values  are  very  small  about  x''=0,  the  normalized  bulk  mean  temperature  on  the  centerline 
is  a  good  way  to  see  the  change  at  the  bulk  temperature  relative  to  the  centerline 

temperature  at  x"^=0,  where  the  bulk  temperature  is  ^  ^  and  x""  is  x""  =  . 

\tg-tg)  RePr 

The  bulk  mean  temperature  vs.  x""  graphs  for  the  Peclet  numbers  1,3,  and  5,  and  the 
thermal  conductivity  ratios  of  1,  1/10  and  the  radius  ratios  of  1.2  and  2  are  seen  in  figure 


21  where  the  bulk  mean  temperature  is 


(te-tp)' 


Figure  20.a  The  normalized  bulk  mean  temperature  with  and  without  wall 
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THE  NORMALIZED  BULK  MEAN  TEMPERATURE 


xPe 


Figure  20.b  The  normalized  energy  absorption  upto  the  entrance  for  parallel  plates  constant  heat  flux  flow 

The  theta  bulk  mean  values  with  the  wall  conduction  case  is  closer  to  the  steady 
state  values  than  the  no  wall  conduction  case  as  it  is  seen  in  figure  21.  The  bulk  mean 
temperature  go  to  zero  as  x  goes  to  plus  infinity.  This  means  that  the  temperature  of  the 
fluid  including  the  wall  conduction  increases  faster  than  the  temperature  of  the  fluid 
neglecting  the  wall  conduction.  The  main  reason  of  the  relative  increase  at  the 
temperature  including  the  wall  conduction  case  is  the  heat  that  is  transferred  by  the  axial 
conduction  of  the  wall  opposite  to  the  flow  direction.  The  heat  moves  to  the  opposite  of 
the  flow  direction  and  is  absorbed  by  the  fluid.  So  the  fluid  temperature  increases  much 
more  than  no  wall  case. 


Figure  21.  The  bulk  mean  temperature  with  and  without  wall  conduction 


2.  Nusselt  Number  vs.  X 


The  Nusselt  numbers  vs.  x+  graphs  are  seen  in  figure  22  for  different  Peclet  numbers, 

/  d 

thermal  conductivity  and  wall  thickness  where  is  x''  = - .  The  Nusselt  number 

RePr 

distribution  is  very  different  from  the  no  wall  conduction  case  as  it  is  seen  in  the  figure. 
As  the  wall  gets  thinner  and  the  thermal  conductivity  increases,  the  values  of  the  Nusselt 
number  gets  closer  to  the  case  of  no  wall  case.  h=l/10  and  k=1.2  is  close  to  the  no  wall 
conduction  case,  so  the  values  of  the  Nusselt  number  for  these  curves  are  close  to  the 
values  of  no  wall  case. 


Figure  22.  The  nusselt  number  distribution  including  the  wall  conduction 
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3.  Heat  Absorption  vs.  X 

The  graphs  of  the  heat  absorption,  E*  vs.  Peclet  numbers  for  the  conductivity  ratios  of 
1,1/10  and  radius  ratios  of  1.2  and  2  were  plotted  in  figure  23  where  the  dimensionless  absorbed 

energy  E*  is  £*  _  ^absorbed  .  As  it  is  seen  from  the  graph.  As  the  wall's  conductivity  or 
-h)  pc  pro 

thickness  increases,  the  heat  absorption  before  x=0  significantly  increases.  The  wall 
conduction  case  has  to  be  added  to  the  problem  for  the  microheat  exchanger  problems. 
Because  the  relative  thickness  of  the  wall  is  very  large  as  it  is  compared  to  the  iimer 
diameter  of  the  tube.  The  values  of  h=1.2  or  h=2  are  realistic  for  the  microheat  machines, 
because  the  diameter  of  the  tubes  are  very  small  for  these  heat  exchanger.  The  heat 
absorption  for  kf7ks=l/10  and  h=1.2  case  is  more  than  ten  times  of  the  case  of  no  wall.  If 
the  wall  thickness  would  be  h=2,  then  the  heat  absorption  would  be  much  greater. 


Figure  23.  The  energy  absorption  upto  the  entrance  with  and  without  the  axial  wall  conduction 
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VI.  CONCLUSIONS 


Axial  conduction  effects  in  duct  flows  are  important  and  have  to  be  considered  for 
cases  of  Pe<10.  The  theta  bulk  mean  temperature  and  Nusselt  number  variations  vs. 
downstream  distance  were  plotted  as  seen  in  the  results.  It  is  obviously  seen  the  increase 
in  the  temperature  against  the  flow  direction  from  x=0  as  the  Peclet  number  decreases. 
Heat  absorption  by  the  fluid  in  the  insulated  region  was  also  quantified  by  using  the  non- 
dimensional  parameter  E*.  The  heat  absorption  increases  logarithmically  as  the  Peclet 
number  decreases,  as  seen  in  the  heat  absorption  graphs.  When  the  wall  conduction 
effects  are  included  in  the  problem,  it  was  calculated  that  most  of  the  heat  absorption  by 
the  fluid  is  because  of  the  wall  conductivity  as  the  wall  conductivity  gets  larger.  The  same 
effect  is  also  seen  with  increasing  wall  thickness.  The  heat  absorption  by  the  fluid  gets 
larger  when  the  wall  gets  thicker.  It  is  very  important  to  include  the  wall  conduction  in 
the  problem  as  the  wall  gets  thicker  and  the  thermal  conductivity  of  the  wall  gets  larger. 
With  large  wall  conductivity  and  wall  thickness,  the  fluid  temperature  increases 
significantly  before  the  heating  starts.  The  current  interests  in  axial  conduction  effects  are 
because  of  the  advent  of  micro  heat  exchanger  technology  with  microtubes  where  the 
Reynolds  number  is  very  small. 
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APPENDIX  -  A  SAMPLE  PROGRAM 


The  Program  for  Circular  Tube,  Constant  Wall  Temperature  Case 

This  program  solves  the  finite  difference  equations  for  the  circular  tube,  constant 
wall  temperature  case  by  using  the  Gauss  Seidel  Iteration  Method  with  the  overrelaxation 
parameter  w=1.8.  The  program  solves  the  equations  doing  the  following  steps  : 

1 .  The  matrix  is  formed  by  the  program,  the  first  row  represents  the  region  which 
touches  the  tube.  The  last  row  is  the  centerline.  The  first  column  is  the  condition,  where  x 
goes  to  minus  infinity.  The  last  column  represents  the  boundary  condition  as  x  goes  to 
plus  infinity. 

2.  The  boundary  conditions  are  applied  to  the  problem.  The  centerline  boundary 
condition,  the  boundary  conditions  at  i=ro,  the  boundary  conditions  as  x  goes  to  minus 
and  plus  infinity  are  applied  to  the  last  row,  fibrst  row,  first  column  and  the  last  column 
respectively. 

3.  The  iterations  are  made  calculating  the  values  at  all  the  grid  points.  As  it  is  seen 
in  the  program,  the  odd  rows  are  calculated  first,  the  even  rows  later.  The  convergence  of 
the  problem  is  better  in  this  case. 

4.  The  difference  between  the  error  and  the  boundary  conditions  for  the  plus  and 
minus  infinity  is  checked.  If  the  boundary  conditions  and  the  difference  between  the 
iterations  are  reasonable,  the  program  is  stopped  and  the  data  is  processed.  The  max. 
difference  between  the  iterations  is  less  than  lE-6.  The  max  difference  between  the 
values  of  the  horizontal  grid  points  which  represents  minus  and  plus  infinity  is  less  than 
lE-6. 
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%  This  problem  solves  the  problem  for  circular  tube 
%  constant  wall  temperature  case 
% 

%  Programmed  by  Ibrahim  Girgin 
% 


format  long 

clear 

start=l; 

iteration_number=300000 ; 

%  Nx  is  the  column  number  of  the  matrix 

%  Ny  is  the  row  number  of  the  matrix 

%  X  zero  is  the  column  number  which  represents  x=:0 

Nx=2501 

Ny=101 

x_zero=1501 

%  For  Peclet  number  is  1 


Pe=l 


%  For  the  relaxation  parameter  1 . 8 
w=1.8 

save_iteration=100 

%  hy  is  the  vertical  distance  between  the  grids 
%  hx  is  the  horizontal  distance  between  the  grids 

hy=l/ (Ny-l) ; 
hx=hy; 


%  teta  is  the  matrix  which  holds  the  grid  points 
%  initialization  of  the  teta 


teta=ones (Ny,Nx) ; 

%  Apply  the  boundary  condition  for  r=Ro,  x>=0 


teta  ( : ,  x_zero  :Nx)  =zeros  (Ny,  length  (x_zero  :Nx) )  ; 


% 

e.g. 

j=Ny 

>  1 

1 

1 

1 

1 

1 

1 

0 

0 

0 

% 

the  initial 

1 

1 

1 

1 

1 

1 

1 

0 

0 

0 

o, 

■« 

value  of 

1 

1 

1 

1 

1 

1 

1 

0 

0 

0 

% 

theta 

1 

1 

1 

1 

1 

1 

1 

0 

0 

0 

a. 

o 

1 

1 

1 

1 

1 

1 

1 

0 

0 

0 

% 

j=l 

>  1 

1 

1 

1 

1 

1 

1 

0 

0 

0 

% 

% 

a- 

i=l 

i= 

:x_zero 

% 

The  values 

for  this  theta 

.  matrix 

is 

. 

0  0  0  0 
0  0  0  0 
0  0  0  0 
0  0  0  0 
0  0  0  0 
0  0  0 


i=Nx 
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o  < 


%  Nx=14  Ny=6  x_zero=8  hx=0.2  hy=0.2 

%  Pe,  w,  file_name, save  iteration  is  what  is  to  be  given 

%  The  variables  that  are  being  used  inside  the  loop 

varl=Pe^2*hx^2; 
var2=2/hy^2; 
var 3=1/2 /hx; 
var4=l/hy^2 ; 
var5=l/varl-var3 ; 
var6=l/varl+var3 ; 
var7=l/2/hy; 
var8=2/varl+4*var4 ; 
var9=2/varl+var2 ; 
varl0=2*hx; 

%  The  iteration  starts  here.  For  every  value  of  kkk,  one 
%  iteration  is  made 

for  kkk=start : iteration  number 


f or_max_err=teta ; 

for  i=2:Nx-l  %  do  the  iteration  for  all  the  columns  except  the 

first 

%  and  the  last  columns 

for  j=l:2:Ny  %  do  the  iteration  for  all  the  odd  rows 

y=(j-l)*hy;  %  calculate  the  y,  0<=  y  <=1 

if  j==l  %  do  the  iteration  for  radius=0  (for  the  center) 

teta  (Ny,  i)  =  (teta  (Ny,  i+1)  *  (var5)  +teta  (Ny,  i~l)  *  (var6)  +  ... 

4  *  teta  (Ny- 1 ,  i )  *  var4 )  /var  8  ; 

elseif  j==Ny  %  do  the  iteration  for  radius=R 

if  i<x__zero  %  if  i<x_zero  &  for  radius  R,  it  is  insulated 

teta (1, i) =  ( (teta (1, i+1) +teta (1,  i- 
1) ) / (varl) +teta (2,i) *var2) /var9; 

end 

else  %  do  the  iteration  for  all  the  odd  rows  except  r=0 

and  r=R 

teta (Ny- j+1, i) =teta  (Ny- j  +1, i) +w* ( (teta (Ny- j+1, i+1) * (var5  + 
y'^2/varl0)  .  .  .  . 
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+teta  (Ny- j+1,  i“l)  *  (var6 
var7/y) + . . . 

teta (Ny- j+2, i) * (var4  - 


-  y^2/varl0) +teta (Ny- j , i) * (var4 
var7/y) ) /var9-teta (Ny- j+l, i) ) ; 


+ 


end 


end 

^or  j=2:2:  (Ny-1)  %  do  the  iteration  for  all  the  even  rows 

y=(j-l)  *hy; 

teta(Ny-j+l,i)=teta(Ny-j+l,i)+w*(  (teta  (Ny- j+1,  i+1)  *  (var5  + 
y^2/varl0) . . . 

+teta  (Ny-j+l,i-l) * (var6  -  y^2/varl0) +teta (Ny- j , i) * (var4  +  var7/y)  + 
teta(Ny-j+2,i) * (var4  -  var7/y) ) /var9-teta (Ny- j+1, i) ) ; 


end 


end 


End  of  the  iteration 


% 

S' 


Save  the  variables  at  every  "save_iteration'' 
e.g.  if  the  save_iteration  equals  to  20,  then  save 

the  variables  at  when  kkk=20, 40, 60, 80 , . . .  so  on. 


if  fix (kkk/save_iteration) ==kkk/save_iteration; 

max__error=  (max  (max  (abs  (teta-for_max_err) )  ) ) 
clear  for_max_err 

save  forpel  teta  Nx  Ny  x_zero  w  Pe  hx  hy  max_error  kkk 
save_iteration  start  y 


end 

end 
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